A mixed population of competing TASEPs with a shared reservoir of particles 
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We introduce a mean-field theoretical framework to describe multiple totally asymmetric simple 
exclusion processes (TASEPs) with different lattice lengths, entry and exit rates, competing for 
a finite reservoir of particles. We present relations for the partitioning of particles between the 
reservoir and the lattices: these relations allow us to show that competition for particles can have 
non-trivial effects on the phase behavior of individual lattices. For a system with non-identical lat- 
tices, we find that when a subset of lattices undergoes a phase transition from low to high density, 
the entire set of lattice currents becomes independent of total particle number. We generalize our 
approach to systems with a continuous distribution of lattice parameters, for which we demonstrate 
that measurements of the current carried by a single lattice type can be used to extract the en- 
tire distribution of lattice parameters. Our approach applies to populations of TASEPs with any 
distribution of lattice parameters, and could easily be extended beyond the mean-field case. 

PACS numbers: 87.10.Hk, 87.10.Mn, 05.40.-a, 05.70.Ln, 05.60.-k 



I. INTRODUCTION 

The totally asymmetric simple exclusion process 
(TASEP) has become a paradigm in non-equilibrium sta- 
tistical physics, due to its rich phenomenology and wide 
applicability pj[6]. The majority of this work has focused 
on the case of a single TASEP with a fixed entry rate (cor- 
responding to infinite availability of particles). However, 
many physical and biological situations involve multiple 
lattices with different parameters, whose dynamics is cou- 
pled via competition for a finite pool of particles. In this 
paper, we present a simple mean-field theoretical frame- 
work to study such scenarios, which can be used for an 
arbitrary number of lattices and for any distribution of 
lattice parameters. We show that non-trivial behavior, 
including extension of the shock phase to a finite region of 
parameter space, and "buffering" of individual currents 
to changes in the reservoir particle number, can emerge 
from the competition between lattices for particles. We 
further show that in these systems, information on the 
whole system can be extracted from measurements of the 
behaviour of an individual lattice subtype. 

The standard TASEP is a stochastic process which de- 
scribes the collective movement of particles along a one- 
dimensional lattice composed of L sites. On the lattice, a 
particle hops stochastically from site i to i + 1, provided 
that site i + 1 is empty. Particles enter the first site 
(if unoccupied) and leave the last site with fixed rates. 
This model shows three distinct phases controlled by the 
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boundary rates: the low density (LD) phase occurs when 
the entry rate is limiting, the high density (HD) phase 
when then exit rate is limiting, and the maximal current 
(MC) phase when both entry and exit rates are large 
enough that the current is limited only by the hopping 
rate along the lattice. Exact steady-state solutions for 
the density of particles on the lattice and the current are 
known [ZHS] . 

First introduced as a model of protein synthesis jTUl 
ITT] , the TASEP nowadays provides a generic model for 
transport processes, with applications including the pro- 
duction of mRNA [H EES] and protein [H [TS] in biolog- 
ical cells, the motion of motor proteins along cytoskele- 
tal filaments 16 -2f]j, the transport of vesicles in fungal 
hyphae [21], collective insect motion [22] and traffic flow 
problems [23J . Most of these studies have focused on par- 
ticle motion on a single TASEP track (including tracks 
with multiple lanes, e.g. [24]); in reality, however, one of- 
ten has multiple tracks, which may have different lengths, 
hopping and boundary rates, and which compete for a fi- 
nite number of particles. For example, in traffic dynamics 
vehicles may distribute themselves across different roads, 
while in biological cells different mRNA molecules com- 
pete for the cellular protein production machinery |25j . 
Our aim in this paper is to provide a simple, intuitive 
and generic framework to study such problems. 

This work is not the first to consider the effects of a 
finite reservoir of particles, or of competition for parti- 
cles between several TASEPs. Recently, Adams et al [26] 
and Cook and Zia [27J studied a single TASEP with a 
finite reservoir of particles, using Monte Carlo simula- 
tions, mean-field theory and both simple and generalised 
domain-wall theories. Cook et al [28] extended this work, 
using Monte Carlo simulations and domain- wall theory to 
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p — a LD phase a < /3, a < 1/2 J = 7a (1 — a) 
p = l-0 HD phase a > /3, /3 < 1/2 J = 7/3(1-/3) 
p=l/2 MC phase a,/3>l/2 J = 7/4 

TABLE I: Summary of the mean-field solutions for a single 
TASEP with fixed entry and exit rate. These relations be- 
come exact in the L — > 00 limit. Note that here a and ft 
denote the dimensionless parameters as defined in the text. 



describe several TASEPs of the same or different lengths, 
sharing the same particle reservoir. This work revealed 
interesting effects due to the reservoir-induced coupling 
between TASEPs. In particular, a regime can exist in 
which reservoir density is independent of the total par- 
ticle number. The domain-wall approach is powerful be- 
cause it can incorporate fluctuations and thus describe 
accurately the density profiles on the lattices. However, 
this approach is difficult to generalize to many lattices 
with arbitrary distributions of entry and exit rates. Here, 
we take a much simpler, mean-field approach. Our ap- 
proach does not provide information on fluctuations and 
cannot predict density profiles. However, it yields a sim- 
ple way to calculate phase diagrams, and is easy to gen- 
eralize to arbitrary populations of lattices. 

In section |TTJ we describe our mean- field approach and 
present relations for the partitioning of particles between 
the reservoir and lattices. We apply this methodology to 
the case of a single TASEP coupled to a finite reservoir in 
section |III A| and to a population of identical lattices in 
section [ill B As a result of the competition for particles, 
a new region of the phase diagram emerges, in which 
the LD and HD phases coexist; the size of this region 
increases as the particle availability decreases. Next, in 
section 



IV we consider a system with two different types 



of lattices and we find that an LD-HD phase transition 
on one lattice type causes the behaviour of the other lat- 
tice type to become independent of the particle number. 
In section [Vj we generalize to the case of mixed popula- 
tions of lattices with arbitrary distributions of boundary 
rates. Here we find that an LD-HD phase transition on 
any lattice type "buffers" all other lattices to changes in 
the particle number, and we show that the entire distri- 
bution of lattice parameters can be extracted from the 
current on a single lattice subtype. Finally we present 
our conclusions in section fVll 



II. 



MEAN-FIELD THEORY FOR MULTIPLE 
COMPETING TASEPs 



We begin our discussion by recalling the behavior of 
a single TASEP with fixed hopping rate 7, entry rate a 
and exit rate /3 [TH5] . For the purposes of this paper we 
discuss only the mean-field results for the steady-state 
average particle density (number of particles per unit 
length) and current on the lattice, which become exact 
in the limit of a large lattice. We first define rescaled 
dimensionless entry and exit parameters a :— a/7 and 



P := (3/-f (i.e. we use the hopping rate 7 to define the 
unit of time) . The (mean-field) density p of particles on 
the lattice takes one of three values, depending on a and 
/3, as summarized in Table [I] The current J is related to 
the density by J = 7p(l — p). In the low density (LD) 
phase, where a < j3 and a < 1/2, the behavior of the 
system is dictated by the entry rate and p = a. In the 
high density (HD) phase, where a > (3 and f3 < 1/2, the 
exit rate is limiting, particles queue on the lattice, and 
p = 1 — j3. When both entry and exit rates are large, 
q,/3 ^ 1/2, neither rate is limiting and the system is 
in the maximal current (MC) phase, for which p = 1/2. 
The transition from either LD or HD to the MC phase is 
continuous in the density, but the transition between the 
LD and HD phases is discontinuous in the density. On 
the HD-LD transition line, where a = f3 and a,/3 < 1/2, 
there is "coexistence" between the LD and HD phases: 
part of the lattice takes on the density of the LD phase 
while the remaining part takes on the density of the HD 
phase, with a well-defined boundary between these do- 
mains, called shock. This scenario is known as a shock 
phase (SP) 0. 




FIG. 1: Schematic illustration of the multi-track TASEP with 
a finite pool of particles. Different lattices may have different 
lengths, entry and exit rates. The entry rate of particles onto 
a given lattice depends linearly on the concentration N r /V 
of free particles, where N r is the number of particles in the 
reservoir, and V is the volume of the reservoir (shaded region) . 



In this paper, we consider the case shown in Figure 
[l] in which multiple lattices with different lengths, en- 
try and exit parameters compete for a finite number of 
particles. The reservoir of free particles is assumed to 
have a finite volume V When a particle exits a lat- 
tice, it immediately enters the reservoir, from where it is 
free to enter any other (or the same) lattice. We assume 
that free particles are well-mixed (not correlated) and 
homogeneously distributed within the reservoir. Since 
particles are conserved, the total number of particles N, 
can be written as 

N = N r + Ni , (1) 
where N r denotes the number of particles in the reservoir 
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and Ni the number of particles on the lattices. In turn, 
Ni can be found by summing over the lattices: 

M 

where M is the total number of lattices, pj is the density 
of particles on lattice j and Lj is the length of lattice j. 

The entry rate of particles onto any given lattice de- 
pends on the density N r /V of particles in the reservoir. 
Since the reservoir is well-mixed, we expect this rate to 
be linear in the particle density, so that an individual 
TASEP j experiences an injection rate aj given by 

N r 

a i :=a o,j-^, (3) 

where ao,j is the "intrinsic" affinity of lattice j for par- 
ticles (with units of volume/time). This relation implies 
that the rate of particle binding events per lattice is pro- 
portional to the free particle concentration: this is ex- 
pected to be true as long as the reservoir does not become 
too crowded with particles. 

In the limit of a large total number N of particles (N — > 
oo at constant N/V), N r = N - N t > N t , so that N r /V 
will tend to the total particle concentration N/V, and we 
recover the standard TASEP with constant entry and exit 
rates. Since the maximum value of TV; is Lj = LM 
(denoting L := (Lj) as the average over all lattices), the 
standard TASEP limit is approached when N ^> LM. 
This means that the finite reservoir plays an important 
role only for N/LM ~ 1 or smaller. Following previous 
work 26-28 , the exit rate from lattice j is assumed to 
be independent of the reservoir, being simply given by 
its "intrinsic" exit rate Pj. 

Within the mixed population of TASEPs, the phase 
behavior of a given lattice j is completely determined by 
its (iV r -dependent) effective entry rate ay and exit rate 
j3j, according to the rules for a single TASEP, as listed 
m Table [TJ For any given value of the reservoir particle 
number N r (assuming fixed V), subpopulations of lat- 
tices will be in the LD, HD and MC phases, depending 
on their values of ay and j3j [33] . We can therefore write 
the total number of particles on the lattices as 

N t = N LD + N HD + N MC , (4) 

where Nld, N^d an d Nmc are the total numbers of 
particles on lattices in the LD, HD and MC phases, re- 
spectively. The density of particles on each lattice de- 
pends on which phase it is in, according to the relations 
m Table [T] Equation @ can therefore be rewritten as 

Ni(N r ) = Yl <*j( N r) L i + E 0- - h) L i + E V2 ■ 

j,LD j,HD j,MC 

(5) 

Here, we explicitly note that Nld depends on the reser- 
voir particle number N r through the dependence of the 
effective injection rate ay on N r . The three sums in 



Eq. (§ are over lattices in the LD, HD and MC phases: 
because ay depends on N r , changes in the reservoir par- 
ticle number will drive phase transitions on the lattices, 
so that these three sums will encompass different lattice 
subpopulations for different values of N r . We can then 
combine Eqs. ([TJ and ^ to express the relation between 
the reservoir particle number N r and the total particle 
number N: 

N = N r +J2 ^{Nr)Lj + E(!- k) L i + E V 2 ■ 

j,LD j.HD j,MC 

Equation ([6]) forms the basis of our theoretical approach; 
in the remainder of the paper we explore its implications, 
first for some simple cases and then for more complex 
scenarios. 



III. A HOMOGENOUS POPULATION OF 
TASEPs COUPLED TO A FINITE RESERVOIR 

We first consider the case of a homogeneous population 
of lattices. Our system thus contains M identical lattices 
of length L, intrinsic injection rate «o and exit rate j3, 
coupled to a finite pool of N particles. Since the lattices 
are identical, in our mean-field approach each lattice ex- 
periences the same effective injection rate a — aoN r /V. 
This implies that all the lattices are in the same phase; 
we refer to this as the "global phase" of the system. The 
relation between N and N r is now given by: 

{N r + LMa N r /V a N r /V < f3, a N r /V < 1/2 
N r + LM(1 - (3) a N r /V > f3, (3 < 1/2 
N r + LM/2 a N r /V ^ 1/2, p ^ 1/2 , 

where the first, second and third relations refer to the 
global LD, HD and MC phases, respectively. 



A. Single TASEP 

For the purposes of illustration, we first consider the 
case of a single (M = 1) TASEP of length L. For this test 
case, our approach corresponds closely to the calculation 
presented for a saturating function a(N r ) by Adams et 
al |26j . The relation between the total particle number 
and the number of particles in the reservoir is given by 
Eqs. 0, taking M = 1. 

Figure^ (solid and dashed lines) illustrates the math- 
ematical mapping between N r and N given by Eqs. ([7|. 
If (3 < 1/2, then when a = (3 the lattice enters the HD 
phase and the number of particles on the lattice jumps 
discontinuously. The results of the mean-field theory 
are in very good agreement with continuous time kinetic 
Monte Carlo simulations for a single TASEP coupled to 
a finite reservoir of particles, for the same parameter sets 
(circles and squares). 
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FIG. 2: (a) Relation between total particle number N and 
reservoir particle number N r for a single TASEP with a finite 
reservoir of particles, from Eqs. Q with M = 1, L = 1000 
and ao/V = 0.001. The solid and dashed lines show the 
mean-field mapping, Eqs. |7| for /3 = 0.4 and 0.6 respectively, 
while the circles and squares show the results of kinetic Monte 
Carlo simulations for the same parameter sets, (b) The same 
data as in (a) but inverted, i.e., showing the reservoir particle 
number 7V r as a function of the total particle number N. 



In most realistic scenarios, however, we expect that the 
total number of particles N is fixed, while the reservoir 
number N r achieves a self-adjusting steady-state value 
depending on N and the parameters of the lattices. We 
therefore need to invert the mapping of Eqs. ^ to find 
N r (N): the number of particles in the reservoir for a fixed 
total particle number. For a single lattice, this inversion 
is shown in Fig. ^p. For f3 < 1/2 the LD-HD phase tran- 
sition occurs over a finite range of N; i.e. there exists a 
range of values of N for which the number of particles 
in the reservoir N r is independent of N. In this range of 
values of N, the lattice is in the shock phase, with HD 
and LD domains. If a new particle is added to the system 
it is quickly absorbed by the lattice rather than staying 
in the reservoir, increasing the size of the HD domain. 
Likewise, if a particle is removed from the system, the 
relative size of the HD and LD domains will adjust to 
keep the number of particles in the reservoir constant. 
Thus, while the lattice is in the shock phase, N r is inde- 
pendent of N: the shock phase lattice buffers the particle 
reservoir to fluctuations in the total particle number. It 
is interesting to note that our mean-field theory captures 
this key qualitative feature of the system, even though 
it deals only with the average density of particles on the 
lattice and cannot resolve the dynamics of the LD-HD 
"domain wall" . 



dent of N. The width of this plateau in N r (N) is exactly 
the height of the step in N(N r ) and can be easily ob- 
tained from Eqs. Q: substituting the critical condition 
iV r = fiV I evo and taking the difference between line 2 and 
line 1 in Eqs. ([7j gives the plateau width 



AN = LM(1 - 2(3) 



(8) 



Hence, AiV increases linearly with the number of lattices 
M. 

We can also use Eqs. ([7| to obtain a global phase dia- 
gram for the system in the parameter space of the intrin- 
sic entry and exit rates (ao/V, /3), for fixed N, L and M, 
by substituting the critical condition N r — jiV/ao that 
defines the phase boundary into the relations for N(N r ) 
in Eqs. ^ [33] . Figure [3] shows this phase diagram, for 
several different values of the total number of particles 
N. Interestingly, as N decreases (Figures [3^i-d) , a finite 
region of the phase diagram emerges where the lattices 
are in the shock phase (SP): here, the particle number is 
too high for the LD phase but insufficient to allow the 
TASEPs to reach the HD phase. As discussed above, in 
this SP region of the phase diagram, the reservoir parti- 
cle number N r remains constant such that aoN r /V — ft. 
The extent of this SP region decreases as the number 
of particles in the reservoir increases, and for N 3> LM 
the phase diagram tends to that of a single TASEP with 
infinite resources and a = aoN/V. 
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B. Multiple TASEPs 

Building on our results for a single TASEP, we next 
consider a system consisting of M identical lattices. This 
system shows qualitatively similar behavior to that illus- 
trated in Fig. [2] for the single TASEP case: if /3 < 1/2, 
the system undergoes a global LD-HD transition when 
aoN r /V = (3, and this produces a range of values of N 
for which the particle reservoir number N r is indepen- 



FIG. 3: Phase diagrams in the (ao/V,j3) plane of the multi- 
track TASEP for a fixed set of parameters (N,L,M). (a) 
N = 1.5 ■ 10 4 , (b) TV = LM = 10 4 , (c) N = 9.5 • 10 3 , (d) 
N = 8 ■ 10 3 . In all cases L = 10 3 and M = 10. These 
choices of parameters represent, as discussed in Sec. |ll] high 
(a), medium (b-c) and low (d) levels of availability of free 
particles. Note that in this figure we consider only cases in 
which 2N > LM. If instead 2/V < LM, there are not enough 
particles to bring all lattices into the MC phase and the phase 
diagram contains only the LD and SP phases. 
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While Fig. [3] provides insight into the physics of the 
system, in a real experimental situation we expect that 
the control parameters are most likely to be the number 
of particles N and the number of lattices M, rather than 
ao and (3. It is therefore useful to plot phase diagrams 
in the parameter space (N,M), for fixed ao/V and (3. 
This can be achieved by reformatting the boundary 
conditions for the global LD, HD and MC phases in 
terms of N and M, using the relations listed in Eqs. tffh. 
For example, for the system to be in the global LD phase 
we require a N r /V < (3 and a N r /V < 1/2. Using 
the relation N = N r + LMa N r /V which holds for 
the LD phase, the first inequality can be reformulated 
as N < (V/3/a.Q + LM/3) while the second becomes 
N < (V/a a + LM)/2. This procedure allows us to 
build phase diagrams in the (N, M) parameter space. 
Figure |4^i shows a case where (3 < 1/2, so that the 
phase diagram contains a shock phase region, while 
Figure [4Jd shows a case where f$ ^ 1/2, so that the 
system instead makes a continuous transition between 
the LD and MC phases. Because of competition among 
the lattices for particles, the global state of the system 
depends on both N and M. Increasing the number 
of lattices M decreases the number of particles per 
lattice, pushing the system towards the LD phase, 
while increasing the number of particles N increases the 
particles per lattice, pushing the system towards the HD 
or MC phase. Table [Tl] in Appendix [7] summarizes the 
boundaries among the global phases of this multi-track 
TASEP, in both parameter spaces (N, M) and (a /V, f3). 
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FIG. 4: Multi-track TASEP phase diagram in the parameter 
space (N, M), for a /V = 5 • 1(T 3 and L = 10 4 . Panel (a) 
shows a case with /3 = 0.25, while in (b) we show a case in 
which the exit rate is not limiting, j3 1/2 (note that in this 
case the boundary between LD and MC does not depend on 
the precise value of 0). 



IV. A MIXED POPULATION OF TASEPs 

We now move on to the more relevant case, where the 
lattices are not all identical. Here, we find that the cou- 
pling between lattices induced by the finite reservoir has 
interesting and non-trivial effects. We first consider a 



population composed of two different types of lattices: 
AfW lattices with intrinsic injection rate and 
lattices with intrinsic injection rate > . Note 
that here we introduce a new notation: we use upper in- 
dexes in brackets (e.g., a^) to indicate properties shared 
by all lattices in the same subpopulation, as opposed to 
lower indexes (e.g., cxqj) which we used to denote the 
properties of individual lattices. For the sake of simplic- 
ity we suppose that all the lattices have the same length 
L and exit rate (3. Our methodology can easily be ex- 
tended to the case of different L and (3 (see Section [vj . 

Because the two lattice subpopulations have different 
values of ao, they will undergo the LD-HD or LD-MC 
phase transition at different values of the reservoir par- 
ticle number N r . For the same value of N r , the two 
subpopulations of lattices can therefore be in different 
phases. The total number of particles can be expressed 
in terms of the particle densities p^ and p^ on the 
type-1 and type-2 lattices: 



N = N r + p {1) LM {1) + p {2) LM {2) , 



(9) 



The densities p^ (i = 1,2) depend on the phase of the 
lattice: if (3 < 1/2 one has 



„(») 



,(0 



if a^Nr/V < 13 
if a^N r /V > /3 
while if /3 ^ 1/2 the densities are 
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(11) 



We can express Eq. ^ in a compact way by making use 
of the Hcavisidc function (defined as 0(z) — 1 for z ^ 
and 6(z) = otherwise). For f3 < 1/2 this results in: 
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which can be rearranged to give: 

N = N r + °„ LM {V) + °„ LM {2) 
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(13) 
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For the case where j3 > 1/2 we could write an equivalent 
equation, using instead the densities and boundary con- 
ditions appropriate for the LD to MC phase transition. 
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It turns out however that this is exactly equivalent to 



Eq. (13 1, but with j3 replaced by 1/2. Equation (13) can 



therefore be used to describe the full behavior of the sys- 
tem, for any value of j3, with the proviso that for f3 > 1/2, 
we simply set j3 — 1/2 in the equation. 
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FIG. 5: Reservoir particle number N r as a function of to- 
tal particle number N for the system with two subpopu- 
lations of lattices (color online), with /V = 5 • 10~ 4 , 
ct^/V = 9 ■ 10~ 4 , (3 = 0.15, L = 300, M (1) = 30, and 
M( 2 ) — 25. The full red line shows the results of the mean- 



field theory, produced by numerical inversion of Eq. ( 13 1; the 



blue circles are the results of kinetic Monte Carlo simula- 
tions with the same parameter set. For the mean-field theory, 
the first derivative of N(N r ) is discontinuous at the critical 
points given by Eqs. (221. The width AiV of the first plateau 
is given by AiV = LM^' 2> (1 — 2/3) while the second plateau has 
AiV = LM (1 '(1 — 2/3), see Eq. Q. The phases of the two lat- 
tices subpopulations are also labelled: for example, LD/HD 
denotes a regime with the type-1 lattices in the LD phase and 
the type-2 lattices in the HD phase. 



From Equation ( 13 1 we can obtain the physical relevant 



relation N r (N) by inversion (as in Sec. Ill B). N r (N) is 
plotted in Figure 5, for the case j3 < 1/2: the two #-steps 
of Equation ( 13 1 appear as two plateaus, corresponding 



to two regimes where the reservoir is buffered due to the 
emergence of a shock phase (SP) on one lattice subtype. 
In the case (3 ^ 1/2 we do not obtain this effect, because 
the particle density on a lattice is continuous across the 
LD to MC transition. Figure [5] also shows the results of 
kinetic Monte Carlo simulations for a multi-TASEP sys- 
tem coupled to a finite reservoir, with the same parame- 
ter set. The agreement between the mean-field solution 
and the simulations is very good; the slight discrepancy 
around the phase transitions is due to finite size effects 
(the mean-field solution becomes exact only in the limit 
L -t oo). 



Equation ( 13 ) also provides a simple way to deter- 



mine the phase boundaries for this system; these oc- 
cur as the discontinuities of the Heaviside step func- 
tions are approached from above and below -i.e. at 



N, 



{v^/afA and N r -> {v(3/a ( ^ . Appendix 
[2] gives explicit forms for these phase boundaries. 

Figure [6] shows that the current per lattice and the 



particle density on the individual lattices exhibit a re- 
markable dependence on the total particle number N. 
For the range of values of N over which lattice type 2 
is in the SP, the current and density on the lattices of 
type 1 remain constant, even though these lattices are far 
away from a phase transition. This buffering effect occurs 
because the reservoir particle number remains constant 
while any lattice subtype is in the SP; this fixes a, and 
hence the current and density, of all lattices which are 
in the LD-phase (see Table 1). Hence, lattices of type 
1 are affected by the phase transition on lattice type 2, 
through the coupling to the finite particle reservoir. Note 
that while lattice type 2 is in the SP, its current remains 
constant, but its density increases linearly with N. This 
reflects the fact that, in the "buffering regime", the SP 
lattices absorb particles as N increases, keeping N r con- 
stant. The same observation applies to lattice type 1 in 
the second buffering region. 
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FIG. 6: Current per lattice J w (a) and density p w (b) on 
type-1 and type-2 lattices for the system with two lattice sub- 
populations. The red and black lines represent the mean- 
field theory results for the type-1 and type-2 subpopulations 
of lattices respectively (color online); the squares and circles 
show kinetic Monte Carlo simulation results. The parameters 



used were a { L> /V = 4.5 • 10" 6 , a ( Q 2) /V = 7 • 10" 6 , /3 = 0.15, 
L = 1000, M (1) = 15, and M (2) = 20. Note that the units of 
current are the particle hopping rate 7. 
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FIG. 7: Multi-track TASEP phase diagram in the parameter 
space (M (1) , M (2) ). Panel (a) presents a case with f3 = 0.25, 
while in (b) we show the situation in which the depletion 
rate is not limiting, /3 ^ 1/2. The other parameters used are 
ai/V = 6.5 ■ nr 3 , a 2 /V = 8.5 ■ 10~ 3 , L = 10 3 , N = 10 4 . 

To conclude this section, we investigate the phase be- 
havior of this system as a function of the numbers 
and M ( 2 ) of lattices in the two subpopulations. The 
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parameter space (AfW, M^) is likely to be relevant 
experimentally: for example, in intracellular transport 
problems, and might represent the number of 

cytoskelctal filaments of two different types along which 
motor proteins can travel. Figure [7] shows the phase di- 
agram of the system in the (M^ 1 ^ M^) plane for (a) 
< 1/2 and (b) ^ 1/2. It is clear that the compe- 
tition for particles plays an important role: the phase 
behavior of each lattice subpopulation depends strongly 
on the size of the other subpopulation. 



V. A MIXED POPULATION WITH 
ARBITRARY DISTRIBUTION OF BOUNDARY 
RATES 



where (L) is the average lattice length and we have in- 
serted the limits of integration as detailed in Table|l] Note 
that since P(j3 > 1/2) = we are free to choose the up- 
per limit of the /3-integration as infinity. This simplifies 
our further calculations. 

After adding and subtracting the term 
(L)M J °° d0 fyp /Nr da P{a , 0) and inserting 
the expressions for the subpopulations into the particle 
conservation equation, we arrive at the result 



N(N r ) = N r + (L)M 



N r (a ) 
V 



(16) 



VP/N r 



We now formulate the concepts of Sections [TI]|IV| in 
a general way, to allow us to describe a mixed lattice 
population with an arbitrary distribution of parameter 
values. As in the preceding discussion - see Eqs. ([!]) 
and ^ - our strategy is to express the total number of 
particles N = N r + Nld + Nhd + Nmc as a function of 
the boundary parameters of the lattices. 

The total number of lattices is denoted by M and the 
normalized distribution of boundary rates and lengths 
(i.e. their relative frequencies) for the mixed lattice pop- 
ulation is denoted by P(«o> Pi L). Such a continuous dis- 
tribution could be appropriate in situations where lat- 
tice parameter values are sensitive to small environmen- 
tal changes, or where the total number of lattices is very 
large. We can then express the total number of particles 
on lattices in the LD, HD and MC phases as: 



N, 



LD 



N 



HD 



mJ dL J J l'^- P{a o ,0,L)da o d0 

LD 

M f dL [[ L (1 - 0) P(ao, 0, L) da d/3 



N, 



MC 



HD 



Mj dL J J -P(a , P,L)dao d/3 , (14) 

MC 



where the integration limits are determined by the con- 
straints on a ~ a N r /V and (3 for the LD, HD and MC 
regions, respectively, as defined in Table [T] Note that the 
MC case can be recovered as a special case of the HD 
with j3 = 1/2. Therefore, we can limit our analysis to 
< 1/2, i.e. assuming P(0 > 1/2) = [35]. 

We also assume that the lengths of the lattices are in- 
dependent of the boundary rates, such that P(ao, Pi L) = 
P(ao, ft)P(L). This allows us to perform the integration 
over L in Eq. (14), yielding 



f°° fNr (YnN 

N LD = (L)M / d0 / da o ^P(a o ,0) 
Jq Jo v 

/>oo />oo 

N HD = (L)M dB da (1 - P) P(a , P) (15) 
Jo J%£ 



where (a<j) = f n °° dB f n °° dap aoP(ao, 0) is the average 
value of otp. Equation ( p~6| ) is a generalisation of Eqs. ^ 
and (13) to a continuous distribution of parameters. 



A. Discrete distributions 

We first consider the case where our mixed popula- 
tion of lattices contains a finite number of distinct sub- 
populations with parameters (ckq , 8^). In this case the 
distribution P(ao,/3) can be expressed as a sum over 8- 
functions: P(a ,B) = £iOT W /M)8(a^ - a a )S(0^ - 
0), where is the number of lattices in subpopulation 
i and M is the total number of lattices. Equation ( 16 ) 
then takes the form 



N(N r 



N r + (L) 



AI 



N r (a ) 
V 



(17) 



N r - 



VB (i) ' 



which is an extension of Eq. ( 13 1 to an arbitrary num- 



ber of subpopulations. The function N(N r ) has discon- 
tinuities at N r = VpM/o$, for B (i) < 1/2. For the 
inverse relation N r (N), these become plateaus of width 
AN = (L)M® (1 - 2B (l) ). Each lattice type for which 
fiw < 1/2 gives rise to a distinct plateau in N r (N). 
Within plateau region i, lattices of type i are in the shock 
phase, at the LD-HD phase boundary. Importantly, be- 
cause the reservoir particle number N r controls the be- 
havior of the whole system, a plateau in N r (N) implies 
that the entry of any subpopulation of lattices into the 
shock phase is sufficient to make the whole system inde- 
pendent of the total particle number N. 

Since the relation N(N r ) cannot be inverted analyti- 
cally, one cannot give a simple prescription for the ranges 
of values of the total particle number where the system 
is buffered. However the generic prescription for the re- 
gions of N over which the system is independent of ./V 



lower boundary of region i 
upper boundary of region i 



lim 

N r ^(Vp(')/a ( l) )- 

lim 



N(N r ) (18) 
N(N r ) 



for all lattice subpopulations i for which /?W < 1/2. 
These boundaries occur at the positions of the steps in 
the 9- functions in Eq. (17); the upper boundaries occur 



as the steps are approached from above, and the lower 
boundaries as the steps are approached from below. It 
is important to note that the phase transitions on any 
given lattice type depend on the parameter values of all 
the lattices in the system, since N r is determined by the 
competition for particles among all the lattices. 



B. Continuous distribution 

We next consider a scenario where the population of 
lattices does not contain distinct lattice subtypes, but 
instead is described by a continuous probability distri- 
bution of lattice parameters P(uo, (3, L). In this case, 
the conservation equation ( 16 ) for the particle number 
does not reduce to a sum of ^-functions, as it does in the 
discrete case, and consequently there are no discontinu- 
ities in the function N(N r ). As a simple example, let 
us assume that all the lattices have the same fixed value 
ao = q:q? with a continuous probability distribution for 
p: i.e. P(ao, (3) = 5(a — o:q)P(/3) (as before we have re- 
defined P(/3) to consider the MC phase as a special case 
of the HD phase). The integral over a in Eq. (16 1 then 
reduces to a step function 9{cxq — (3V/N r ), leading to: 



N(N r ) = N r + (L)M 



a*N r 
V 



(19) 



]N r /V 



1-/3- 



a* N r 
V 



where the upper limit on the integral reflects the condi- 
tion a > f3 for the HD phase. 



To explore the consequences of Eq. (19), we consider 



the specific case where the distribution of exit rates is 
Gaussian: P(f3) = 1/V27rer 2 exp(-(/3 - (/3)) 2 /2cr 2 ). In 
this case, the integral in Eq. ( 19 1 can be calculated ana- 
lytically to give: 



N(N r ) 



Nr 



a* N r 
V 



(20) 



(L}M 



1 - 



v 



Figure [8^, shows the inverse relation N r (N), for two dif- 
ferent widths a of the distribution of (3 values. The most 
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FIG. 8: (colour online) Particle reservoir concentration as 
a function of the total number of particles for normal dis- 
tributed exit rates (a) and exit rates following a bimodal 
distribution a bimodal distribution (b). In both cases, L = 
300, M = 30, a^/V = 10 -4 . In panel (a) the average value 
of the Gaussian peak is chosen (/3) = 0.2, the dashed (black) 
line represents the distribution with a — 0.05, the full (red) 
line represents the distribution with a — 0.01. In (b) the 
two peaks are centered at /3i = 0.15 and {3% = 0.3, respec- 
tively. The dashed (black) line represents the distribution for 
a = 0.025, the full (red) line for a = 0.005. 



striking feature is the "quasi-plateau" at intermediate 
values of N, which mimics the true plateaus observed 
for the discrete case (see for example Figure [5]). As the 
distribution of (3 values narrows (decreasing a) this fea- 
ture becomes closer to a true plateau. The fact that 
this "quasi-plateau" in N r (N) is observed for a rather 
generic continuous distribution of (3 values suggests that 
the buffering of the particle reservoir by lattices entering 
the shock phase is a general phenomenon, with smooth- 
ing of the distribution of lattice parameters tending to 
"soften" the buffering effect. For large N (which implies 
large N r ), the terms in Eq. (20) which are linear in N r 
dominate the /3-dependent terms so that N r (N) becomes 
linear and independent of the /3-distribution. 

Figure [SJd shows the corresponding results for a bi- 
modal distribution of /? values, consisting of two Gaus- 
sian peaks. Once again, for narrow Gaussian peaks, a 
plateau-like form for N r (N) is recovered, but this time 
with two "quasi-plateaus" . This is analogous to the case 
studied in Section |IV[ where each lattice subpopulation 
produces its own range of values of N over which the 
system is buffered. 



C. Relation between the distribution of boundary 
rates and single lattice properties 

A key feature of the systems discussed in this paper 
is that the behavior of all the lattices is coupled via 
the shared particle reservoir, so that the function N r (N) 
depends on the entire distribution of lattice parameters 
P(cto,f3). This implies that measurements of the reser- 
voir particle number N r , or of the current on a few lat- 
tices, as functions of N, contain information on the full 



9 



distribution of lattice parameters P(ao,(3). In this sec- 
tion, we briefly sketch how such measurements could be 
used to compute P{ao, (3). For simplicity, we assume that 
the entry rate ao = is fixed for all lattices, so that our 
aim is to compute the distribution of the exit rates P(/3) . 

We first consider the case where one is able to mea- 
sure experimentally the function N r (N). In this case, 
P{0) can be obtained from the derivative dN r /dN: dif- 
ferentiating Eq. ( [T9| with respect to N r produces a linear 
differential equation for the cumulative probability distri- 
bution Q{P) := fPp(/3')dP'. This differential equation 



dN 
dN~ r 



(L)Ma* 



(L)M 



V 



V 



(21) 



(l-2/3)Q'(^)-^Q(/3) 



where the dependence of Q on (3 is given by substitut- 
ing a^Nf/V = (3. Note that here we used Leibniz' in- 
tegral rule ^Jof(x,y)dx = f{y,y) + J Q y ^f{x,y)dx, 
with x = j3', y — a* N r /V and the integrand f(x,y) — 
(1 — x — y)P{x). Equation (21 1 depends on dN/dN r — 
1 / (dN r /dN) and can be solved by standard methods (e.g. 
the method of variation of constants). If, rather than 
knowing N r (N), we know the current Jj on a particular 
lattice j (in the LD phase) as a function of N, we can 
use the relation dJj/dN = (dJ j / dN r )(dN r / dN) (since 
Jj depends only on N r for fixed V, M^ l \ and fij) to 
write 



dN 
dW r 



dJ 3 

dN 



and note that dJj/dN r can be computed from the 
TASEP result (for lattices in the LD phase) Jj = ja.j(l — 
ctj) where ctj = a^Nr/V. Having thus obtained dN/dN r , 
we can again use the derivative of Eq. ( [l9| to extract 
P(f3). Note, however, that this procedure only works for 
the range of values of (3 for which lattice j remains in 
the LD phase; if the lattice is in the HD or MC phase, 
the current Jj(N) contains no information on the reser- 
voir particle number and cannot be used to obtain P{j3). 
While it remains to be seen how useful the prescription 
outlined here would actually be for extracting P(f3) from 
real (noisy) data, this discussion highlights the impor- 
tant point that, in principle, one can extract information 
on the parameter distribution of the whole system from 
measurements of the behavior of just a single system com- 
ponent. 



VI. DISCUSSION 

In this paper, we have presented a mean-field theo- 
retical framework to study systems in which multiple 
TASEPs with different parameters compete for a com- 
mon pool of particles. We expect this approach to be 



useful in modelling a wide range of systems, from control 
of gene expression in biological cells to traffic flow prob- 
lems. Previous work has addressed the effects of a finite 
particle reservoir on TASEP dynamics, for single lattices 
[2lfl [2"7] and for multiple lattices with equal boundary 
rates 28 ); here, we extend this work to mixed popula- 
tions of lattices with an arbitrarily complex distribution 
of parameters. 

Our theoretical approach, presented in Section [TTJ, is 
based on combining the equation for conservation of the 
total number of particles with the mean-field results for 
the standard TASEP. This approach provides a simple 
way to deal with mixed populations of lattices. Although 
our mean-field theory does not provide information on 
fluctuations or on density profiles, it nevertheless reveals 
interesting phenomena which emerge from the competi- 
tion for particles in a mixed multi- TASEP system and 
provides a method to calculate the full phase diagram. 
These phenomena arise because the finite reservoir effec- 
tively couples all the lattices, so that a phase transition 
on one lattice influences the behavior of the others. Al- 
though in this work we used the mean-field TASEP re- 
sults, one could easily incorporate into the same frame- 
work exact or simulated relations for the particle density 
as a function of the entry and exit rates. 

A key observation which emerges from our work is that 
for a mixed population of lattices coupled to a finite reser- 
voir of particles, any lattice subtype which enters the SP 
absorbs all further particles added to the system, buffer- 
ing the particle reservoir and making the currents and 
densities of all other lattices insensitive to changes in the 
total particle number N, even though these lattices may 
be far from their phase boundaries. This effect is spe- 
cific to systems where the lattices have different intrinsic 
entry and exit rates: if Qo and j3 are the same for all 
lattices (and if lattices are large enough to neglect finite 
size effects) , phase transitions occur on all lattices at the 
same critical particle reservoir number and the current 
does not show an additional plateau as a function of N. 
This was the case in previous work |26j , which considered 
a mixture of lattice with different lengths but identical 
boundary rates. In these conditions, a plateau in the cur- 
rent as a function of N is due only to finite size effects. 

The physical mechanism underlying the buffering of 
the system to changes in N is attributable to lattices un- 
dergoing the LD-HD transition that "soak up" changes 
in the reservoir particle number N r . As they undergo 
this transition, lattices enter the shock phase, in which 
a queue of particles forms at the end of the lattice. For 
shock phase lattices, the particle entry rate a = aoN r /V 
is fixed by the exit rate (3 [a = 0): thus the number of 
particles on a shock phase lattice adjusts to compensate 
for changes in N r . In the phase diagram for the stan- 
dard TASEP model, the shock phase occurs only on the 
line separating the HD and LD phases, where a = /3; 
in the case of a finite reservoir of particles, however, the 
shock phase occupies a finite region of the (ao,/3) phase 
diagram. This is because the same entry rate a can be 
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achieved over a range of N, N r being set by the position 
of the domain wall. 

An interesting analogy can be drawn between this phe- 
nomenon and first order phase separations in (equilib- 
rium) thermodynamics. The plateau in N r (N) which 
arises in our models is a direct consequence of the discon- 
tinuity in the particle density as a lattice undergoes the 
LD-HD phase transition. Similarly, a first order phase 
transition such as the boiling of water involves a disconti- 
nuity in the entropy, which is associated with latent heat: 
during the transition, the temperature remains constant 
even though further heat energy is constantly being sup- 
plied. In this analogy, heat plays the role of N in our 
models while temperature plays the role of the reservoir 
particle number N r . 

We also show in this paper that the coupling between 
lattices induced by a finite particle reservoir makes it pos- 
sible (under some circumstances) to extract the entire 
distribution of lattice parameters from measurements of 
the reservoir density, or indeed the current carried by a 
single lattice subtype, as a function of the total particle 
number N . It will be interesting to explore the feasibil- 
ity of this approach for extracting information from real, 
noisy, experimental data. 

The multi-track TASEP with finite particle reservoir 
studied here bears a similarity to previous work on 
TASEPs on closed networks 29J. In particular, a system 
of M lattices with a common reservoir can be mapped 
onto a network topology formed by M rings having a 
unique common site |30j . The common site plays an 
analogous role to a particle reservoir. However, in that 
problem, in contrast to the one studied here, both the 
exit and entry rates depend on the occupation of the 
common site. 

A key priority for future work must be to explore ways 
to include the effects of fluctuations, which arc neglected 
in our mean-field approach. Previous work has shown 
that interesting fluctuation-driven effects, including lo- 
calisation of domain boundaries, can occur in TASEPs 
with finite particle number [27]; extending this work to 
complex mixtures of TASEP is likely to prove fruitful. 
Another promising avenue may be to study how the be- 
havior of systems of the type studied here changes with 
changes in their volume: this should prove relevant when 
modelling transport or protein production dynamics in 
growing biological cells. 

In summary, we have presented a simple and intu- 
itive mean-field theoretical framework for studying multi- 
TASEP problems with finite reservoir of particles. Our 
method has allowed us to show that interesting physi- 
cal phenomena emerge from the competition for parti- 
cles among non-identical lattices, including buffering of 
the system to changes in the total particle number. This 
approach should prove a versatile tool for studying a wide 
variety of "real-world" problems [31] involving competi- 
tion among complex populations of transport processes. 



Acknowledgments 

We thank Chris A. Brackley, Michael E. Cates, Martin 
R. Evans, Marco Thiel, Ian Stansfield and the anonymous 
referee for valuable discussions and comments on the 
manuscript. RJA was supported by a Royal Society Uni- 
versity Research Fellowship, LC by a SULSA studentship 
and partially by the GDRE 224 GREFI-MEFI CNRS- 
INdAM; PG by a DAAD postdoc fellowship and by EP- 
SRC under grant EP/E030173, and MCR by SULSA and 
BBSRC (BB/F00513/X1, BB/G010722). The collabora- 
tion leading to this work was facilitated by the StoMP re- 
search network under BBSRC grant BB/F00379X/1 and 
by the e-Science Institute under theme 14 "Modelling and 
Microbiology" . 



Appendix: Critical points 

In this appendix, we provide explicit expressions for 
the location of the critical lines for the homogeneous and 



mixed populations of TASEPs discussed in Sections III B 
andllVl 



1. Homogeneous multi-track TASEP 

We first present the phase boundaries for the homoge- 
neous case of Section |III B| in which the M lattices have 
identical boundary rates ag and /3. In our mean-field ap- 
proach, all the lattices are in the same phase, which we 
denote the "global phase" of the system. Table [TT] gives 
the conditions determining the global phase diagram, for 
the cases where the exit rate (3 is limiting (/3 < 1/2) and 
where j3 is not limiting (f3 ^ 1/2). As discussed in Section 
|IIIB[ one may choose as variable parameters either the 
entry and exit rates a® an d /?, or the number of particles 
N and lattices M. These alternative choices correspond 
respectively to the phase diagrams shown in Figures [3] 
and |3J Table [TT] gives explicit forms for the phase bound- 
aries in both these cases: the upper line in each row gives 
the conditions on cxq and j3 (assuming fixed N and M), 
while the lower line gives the conditions on N and M (for 
fixed q;o and ft). 



2. Mixed population of TASEPs 

For the mixed multi-track TASEP discussed in Section 
WJ which is composed of two subpopulations of lattices 
with different intrinsic entry rates ao, the states of the 
system are characterised by the phases of the two lattice 
subpopulations (e.g. in the LD/SP state, subpopulation 
1 is in the low density phase while subpopulation 2 is 
in the shock phase) . Boundaries between these states oc- 
cur when the ^-functions in Eq. ( 13 ) are approached from 



above or below (at these points one of the lattice subpop- 
ulations undergoes a phase transition). More precisely, if 
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TABLE II: Explicit expressions for the constraints to be sat- 
isfied in each global phase of the homogeneous multi-track 
TASEP with M lattices of length L, each with intrinsic en- 
try rate ao and exit rate /3, with TV particles. The phase 
boundaries are given in the case (3 < 1/2 (left column) and 
f3 ^ 1/2 (right column). Each row corresponds to a global 
phase (LD;HD;MC;SP) and for each phase the constraints 
are given in terms of ao and ft (upper line; corresponds to 
Figure [3| and N and M (lower line; corresponds to Figure HJ. 



P < 1/2 



P > 1/2 



< 



P 



LD 



V N - LM/3 



B 

N < V — + LM8 



2JV - LM 



V LM 
N < + 

2a a 2 



P 



LM(1 - 8) 



HD 



MC 



V > N - 
N > V— + LM(1 - 8) 



— > 

V IN - 



N > 



2a a 



LM 
LM 
2 



P 



SP 



N - LM8 V N - LM(1 - 8) 

P P 
V— + LMB < N < V — + LM{1 - 8) 



f3 < 1/2 the phase transitions are located at: 



LD/LD - 


LD/SP : 


N r - 


lim 


N(N r ) (22) 


LD/SP - 


> LD/HD : 


JV r — 


lim 

{vp/^y 


N(N r ) 


LD/HD - 


■> SP/HD : 


N r - 


lim 

(vp/^y 


N(N r ) 


SP/HD - 


> HD/HD : 


N r - 


lim 

{vp /a py 


N(N r ) , 



while if ^ 1/2: 

LD/LD LD/MC : 
LD/MC MC/MC : 



lim 

>y/2, 

lim 



N r ^V/2a { 2) 



N r ^V/2a (2) 



N(N r ) 

N(N r ) 



Tables III and IV give explicit expressions for the con- 
straints on the parameters N, L, V, M^ l \ M^ 2 \ 
and /3, for each of the possible states of the system, 



in the cases where /3 < 1/2 (Table III) and /3 > 1/2 



(Table IV I. To obtain these phase boundaries we insert 



the limits defined in Eqs. (22) and p3l into Eq. (13 1, 



using lim z _j. + 9(z) = 1 and lim 2 ^ - 8(z) = 0. Note 
that this simply means to substitute the critical values 

TABLE III: Phase boundaries for the multi-track TASEP 
with two lattice subpopulations introduced in Sec.|IV| for the 



case p < 1/2, obtained by combining Eqs. ||22[) with (| 13 



Low Density/Low Density 

VP "i 11 , n r,i 
N < —^ + P^LM^ +/3LM< 2 » 



Low Density/Shock Phase 



+ P%yLMW + PLMW < JV < ^§ + P^LM' 1 ' + (1 - /3)LM< 2 > 



Low Density/High Density 

^£r+ P^LM (1) + (1 - (j)LM {2) < JV < J T ^ T + /3LM (1) + (1 - 8)LM 

Shock Phase/High Density 



-^jy + PLM {1) + (1 - P)LM {2) < JV < v -jYj+ (1 - P){M {1) + M {2) )L 



High Density /High Density 



vp 



JV > -tjt+ (1 - P){M (1) + M (2) )L 



N r = pV/af (LD-HD transition) and N r = Vj2af 
(MC transitions) respectively, while the limit from below 
corresponds to substituting the ^-function in Eq. ( 13 ) 
with 9(z) = for the lower limit and 9(z) = 1 for the 
upper limit. 



TABLE IV: Phase boundaries for the multi-track TASEP 
with two lattice subpopulations introduced in Sec. |IV| fo r the 
case p ^ 1/2, obtained by combining Eqs. (23 1 with (131. 



Low Density/Low Density 



LM 



12, 



(23) Low Density/Maximal Current 



V a 



LM (2) 



+ -^tLM (1) + — < JV < 



V L(M (1) + M (2) ) 



Maximal Current /Maximal Current 

V L(M< 1 >+M< 2 ') 



JV > 
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